F = @(x,T) -T.*log(1+exp(-x/T));

vx = linspace(-4, 4, 100);
T = 1.0;
ve = 1./(vx.^2 + 1) + 0.5;
ve0 = 0.5*ones(1, length(ve));
veF0 = F(0.5, T)*ones(1, length(ve));
vr = F(ve, T);

close(figure(1)); figure(1);
plot(vx, vr, '-ro');
hold on;
plot(vx, ve, '-bo');
plot(vx, ve0,  '-g');
plot(vx, veF0, '-g');

Num = 1200;
L = 40;
vx = linspace(-L/2, L/2, Num + 1);
vnx = (-Num/2 + 1):(Num/2);
vnp = (-Num/2 + 1):(Num/2);
disp(norm(vx(2:end) - L/Num*vnx))

vp = 2*pi/L*vnp;

mx = ones(Num, 1)*vnx;
mp = transpose(vnp)*ones(1, Num);

U = 1/sqrt(Num)*exp(2*pi*1j/Num*(-mx.*mp));
disp(norm(U'*U - eye(Num)));

vx = vx(2:end);
T = 1.0;
ve = 1./(vx.^2 + 1) + 0.5;
ve0 = 0.5*ones(1, length(ve));
veF0 = F(0.5, T)*ones(1, length(ve));
vr = F(ve, T);

vr0 = vr - sum(vr)/Num;
vr0p = transpose(U*transpose(vr0));
close(figure(2));
figure(2);
plot(vp,vr0p,'-ro');

close(figure(3));
figure(3);
plot(vx, vr0, '-bo');
hold on;
plot(vx, vr, '-go');


Num = 100;
vm = 1:Num;
vn = 1:Num;
mm = ones(Num,1)*vm;
mn = transpose(vn)*ones(1,Num);

M = 2/pi./(mm + mn);
[v,d] = eig(M);
disp(norm(v'*M*v - d));

Num = 1200;
L = 40;
vx = linspace(-L/2, L/2, Num + 1);
vnx = (-Num/2 + 1):(Num/2);
vnp = (-Num/2 + 1):(Num/2);
disp(norm(vx(2:end) - L/Num*vnx))

vp = 2*pi/L*vnp;

mx = ones(Num, 1)*vnx;
mp = transpose(vnp)*ones(1, Num);

vpa = 1:10;

vx = vx(2:end);

T = 1.0;
ve = 1./(vx.^2 + 1) + 0.5;
ve0 = 0.5*ones(1, length(ve));
veF0 = F(0.5, T)*ones(1, length(ve));
vr = F(ve, T);

for ip = 1:length(vpa)
    vpa(ip) = fn_pa(ip, vr, vx);
end

[M, v, vd] = fn_M(length(vpa));

vpa_t = transpose(v'*transpose(vpa));

[M,v,vd] = fn_M(2);






function [M, v, vd] = fn_M(Num)

vm = 1:Num;
vn = 1:Num;
mm = ones(Num,1)*vm;
mn = transpose(vn)*ones(1,Num);

M = 2/pi./(mm + mn);
[v,d] = eig(M);
vd = diag(d);

end


function [M, v, vd] = fn_M2(Num)

vm = 1:Num;
vn = 1:Num;
mm = ones(Num,1)*vm;
mn = transpose(vn)*ones(1,Num);

M = 2/pi./(mm + mn).*(2 + mm + mn)./(mm + 2)./(mn + 2);
[v,d] = eig(M);
vd = diag(d);

end


function res = fn_pa(n, vf, vx)

va = fn_a(n, vx);

v_nt = zeros(1,length(vx));
v_len = vx(2:end) - vx(1:end-1);
v_nt(1:end-1) = v_len;
v_nt(2:end) = v_nt(2:end) + v_len;

res = sum(v_nt.*va.*vf);

end


function vf = fn_a(n, x)

vf = 1/2/pi*n./(x.^2 + (n/2)^2);

end
